Nonequilibrium electron transport in strongly correlated molecular junctions 
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We investigate models of molecular junctions which constitute minimal Hamiltonians to account 
for zero-bias-anomaly and the satellite features of inelastic transport by molecular phonons. Through 
nonlinear transport calculations with the imaginary-time nonequilibrium formalism, a HOMO- 
LUMO model with Anderson-Holstein interaction is shown to produce co-tunneling conductance 
peak in the vicinity of Kondo resonance which is mediated by a re-emergent many-body resonance 
assisted by phonon excitations at bias equal to the phonon frequency. Destruction of the resonance 
leads to negative-differential-resistance in the sequential tunneling regime. 
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Strong correlation in nonequilibrium electron trans- 
port has emerged as one of the most exciting fields of 
condensed matter physics. The research in this field so 
far has been mostly driven by semiconductor-fabricated 
quantum dots (QDs). The zero-bias anomaly (ZBA) phe- 
nomena have been extensively studied in the context of 
Kondo phenomena [H, In recent years, similar ZBA 
phenomenon in molecular junctions 3, 4, 5, 6] has gener- 
ated tremendous excitement for possible different mech- 
anisms for strongly correlated transport. 

Currently, the research on molecular junctions in 
strongly correlated regime is, both experimentally and 
theoretically, at an early stage and little is known for the 
underlying transport mechanisms. One of the most out- 
standing transport phenomena in molecular devices is the 
co-existence of the ZBA and the inelastic conductance 
peaks, presumably due to molecular phonons [1, 0, Hi- 
Theoretically, the strong correlation in molecular sys- 
tems poses a great challenge since strong Coulomb and 
electron-phonon (el-ph) interactions make perturbative 
approaches unreliable. Only recently, strong correlation 
physics in the Anderson-Holstein model has been under- 
stood for equihbrium systems^, H, @|- Most works on 
nonequilibrium transport in molecular systems have been 
perturbative and often excluded Coulomb interaction [lo'. 
Although nonperturbative nonequilibrium theories 
have seen important breakthroughs 12, 1^ 14, 15, in 
the past few years, the methods have not been adequate 
to tackle complex models such as molecular junctions. 

The main goal of this work is to identify minimal 
Anderson-Holstein models which can describe the Kondo 
anomaly and the inelastic features at finite source-drain 
bias, and reproduce some of experimental findings 0, 0, 
S 0]- We apply the recently developed imaginary-time 
theory and numerically solve the Anderson-Holstein 
models via quantum Monte Carlo (QMC) method [l7| . 
Due to the diverse molecular systems, it is very impor- 
tant at this stage to have guiding principles to catego- 
rize molecular models for different transport phenom- 
ena. The main system of focus here are molecular quan- 
tum dots which exhibit the ZBA accompanied by con- 



FIG. 1: Schematic energy diagrams for isolated molecular 
configurations with respect to source/drain reservoirs under 
bias $. (a) Single-orbital model with one electron occupying 
the level at energy e^- Phonon excitation level is marked by 
a dashed line, (b) With an extra electron, the energy level is 
pushed up ed + U ^ U/2 hy Coulomb repulsion, (c) HOMO- 
LUMO model with the level spacing A. (d) Charge-excited 
state. With A > U/2, the LUMO level is nearly empty. 



ductance oscillations at bias near the ZBA energy scale, 
which have been often attributed to the molecular vibra- 
tions. This problem, from the strong-correlation point of 
view, is quite puzzling since near the Kondo anomaly the 
charge fluctuations are strongly suppressed and phonons 
which interact with electric charge fluctuations are ef- 
fectively decoupled 



18[. In single-orbital Anderson- 



Holstein models, it has been shown that phonon spectral 
features in strong Coulomb limit are weak. 

To resolve the issue, we consider two scenarios. First, 
we note that, at finite bias, the strong correlation ef- 
fects become weaker, as shown in the disappearance of 
Kondo peak [H, Then incoherent charge fluctua- 

tions induced by nonequilibrium may enhance the effec- 
tive el-ph interaction. Within this scenario, we study 
the single-orbital (SO) Anderson-Holstein model. [See 
FIG. [ija-b)] Second, we study a two-orbital model with 
highest-occupied-molecular-orbital (HOMO) and lowest- 
unoccupied-molecular-orbital (LUMO). Despite being 
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more realistic, HOMO/LUMO (HL) models have not 
been extensively studied due to their complexity. Multi- 
ple orbitals allow the electron density distortion to couple 
to molecular distortions, i.e. molecular Jahn- Teller (JT) 
modes, without invoking on-site charge fluctuations. It 
has been shown that in strong correlation limit Q, the 
JT couphng becomes very effective. [See FIG. [l{c-d)] 

For both models, the electron source/drain reser- 
voirs are modeled by the Hamiltonian He = 
Sfeacr ^akacl^ka^aka in terms of the electron creation (an- 
nihilation) operator cl^i;„{caka) with k the continuum in- 
dex, a spin index, and the reservoir index a = ±1 for 
source (L) and drain (i?), respectively. The charge part 
of the QD Hamiltonian for the single-orbital model is 



Hch,SO — Ed ^ 



nda + ^iridT +ndi-lf, (1) 



with the number operator nda = d^^dcr of the QD orbital 
dj., the level energy €d and the Coulomb parameter U. 
The Holstein phonon and the el-ph coupling for the SO 
model can be written as 

Hph,so = t^pho^a + gep{a^ + a){nd] + Udi - 1), (2) 

with for creation of phonon, uipi^ the phonon frequency, 
(7ep the el-ph coupling constant. The tunneling part is 
Ht,so = -tHakai^lcaka + h.c.) with the hopping pa- 
rameter t. Here the tunneling rate is parametrized by 
the hybridization function Tl,r = TTt'^N{0) with A^(0) 
the density of state of the reservoirs. Throughout this 
work, we assume Fl = r_R and use r = rL+Fjf = l 
as the unit of energy. The total Hamiltonian is H — 
Hc + Ht + Hch + Hph. Here, we study the regime where 
the phonon frequency is comparable to the Kondo tem- 
perature and we set Uph ~ F, in contrast to semiconduc- 
tor QD models where the phonon-excited QD-levels are 
discrete and well-defined {tOph > F) [19]. 

We solve steady-state nonequilibrium using the 
imaginary-time formalism [l3| . This method combines 
the nonequilibrium quantum statistics and quantum dy- 
namics within the equilibrium theory via an imaginary- 
time Hamiltonian with complex chemical potentials 
parametrized by the Matsubara voltage ipm — A'KmT as 

k{i^^) ^ ko{iv„,) + V = Ho + ]^{iV^-^)% + V, (3) 

where the many-body interaction is given by V and the 
non-interacting part by Hq = + H^ + ed X^o- ''^da- Pop- 
ulation of the scattering states for source and drain in the 
non-interacting limit is imposed by the operator jT^ . [20j 
^0 = Jlkai'^Lka'^Lka- " i^Rka'^Rka), with the Scattering 
state operator V'lfea from the a-reservoir [23|. Yq can 
be exactly solved in the non-interacting limit. In a per- 
turbation expansion with y, the quantum statistics is 
unaffected by icpm due to e~^^° = e~^'^^°~^^°\ since 



^/3(pm^o represents 2tt x (integer) with respect to the un- 
perturbed scattering-state basis. The quantum dynam- 
ics, represented by an energy denominator in Green func- 
tions, is recovered by the analytic continuation iipm — > 
This formalism can be shown to be equivalent to the re- 
tarded Green function in the Keldysh formalism. 

With this formalism, the equilibrium auxiliary-field 
QMC method [l3| can be immediately applied with 
K{i(pm) as the Hamiltonian. The resulting self-energy 
S(ia;„,i(^m) at the fermion Matsubara frequency tOn = 
{2n + 1)'kT should be analytically continued numerically. 
This is achieved by making an ansatz on the spectral 
representation [l3l |. 



(T-y(e)(ie 



(4) 

with the spectral function cr.y(e). The index 7 of odd 
integer is a combination of reservoir indices a in particle- 
hole lines in a self-energy diagram. This representation 
is exact in the equilibrium limit or in the second order 
perturbation in nonequilibrium. We use o'.^(e) as fitting 
parameters to the numerical self-energy. In the particle- 
hole symmetric case, the a;„-independent term a{iipm) is 
zero. With particle-hole asymmetry, we use a simple-pole 
approximation a{iip„i) = ao + ai[{iipm — z)~^ + {—iipm — 
z*)^^], with fitting parameter oq, ai, z with Imz < for 
ip„i > 0- o,{i^m) has a weak dependence on (p.^^ and the 
analytic continuation has been insensitive to the choice 
of a fitting form. Once all the fitting parameters are 
found, we set i(pm $ and iw„ lo + ij], and obtain 
the retarded self-energy and QD Green function G"''^*(ti;). 
The current is calculated from 



/ = ^ / dujirTAiu;) 



, (5) 



with A{lj) = -7r-iG"''=*(w) [21|. 

The differential conductance G — e(dl/d^) for the SO 
model is shown in FIG. in the particle- hole symmet- 
ric limit Cd = —U/2. At U = 0, the system has a ZBA 
peak with the HWHM {^hwhm ^ 0.3), much reduced 
from the non-interacting value ^ %,u/- unf = 2, demon- 
strating the charge-Kondo effect (9l.l22l| . As U approaches 
the el-ph binding energy 2g'^p/uiph, the correlation ef- 
fect becomes weaker indicating the competition of the 
attractive el-ph and repulsive Coulomb interactions. As 
U grows further, the system approaches the usual spin- 
Kondo regime. The fine structure in the conductance 
shows faint oscillations, reminiscent of some of the exper- 
iments 0. However, given the numerical uncertainties, 
we cannot conclude that the features are a direct man- 
ifestation of inelastic excitation of phonon quanta. In 
an extensive set of calculations we found no well-defined 
phonon satellites near the ZBA energy scale. To fur- 
ther support the idea, we doubled the phonon frequency 
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FIG. 2: (Color online) Differential conductance in single- 
orbital Anderson-Holstein model, (a) Conductance with large 
el-ph coupling {2g1j,/Ljph > U) in the charge-Kondo regime. 
As U increases from zero, correlation effects become weaker, 
(b) Spin-Kondo regime recovered with large U. (c) Larger 
phonon frequency at fixed a^rtluiyh ~ 1. (d) Comparison of 
the Qep = limit to Ref [14j and the negative U model. The 
unit of energy is the non-interacting line- broadening F. 



{ujph = 2) at fixed glp/ujph- YlG.\^c) shows similar re- 
sults as (a-b) with weaker fine structures, possibly from 
more efficient QMC sampling at high ujph- We conclude 
that the main effect of el-ph interaction in the SO model 
is the reduction of the Coulomb interaction and that the 
conductance at high bias is more related to the el-ph 
scattering than to the energy exchange with phonon. 

FIG. m^d) shows results from pure Anderson mod- 
els. In the previous work of Han and Heary the 
particle-hole symmetry condition on the spectral function 
in Eq. (|4]), cf~f{e) = (t_,^(— e), was not properly imposed 
and they obtained underestimated ZBA peak-widths. A 
modified fit gives an improved agreement with Ref. [l^ . 
We also test the idea whether the nonequilibrium im- 
posed on the charge variable as opposed to the spin vari- 
able has any significant effects on the nonlinear transport. 
The positive-negative U models are interchangeable in 
equilibrium [23 | by switching the charge and spin vari- 
ables. The calculation shows that the difference between 
the two models even at high bias is minimal. 

We now turn to the HOMO/LUMO model. We denote 
the QD levels by with i = 1, 2 for HOMO and LUMO, 
respectively. The charge part of the Hamiltonian is 



H, 



U 



h,HL = ^ Q"-i<T + y ("IT + nil - 1)^' (6) 



with the HOMO level at ei = e^, the LUMO level at 
£2 = + A. The Coulomb interaction is set to act only 
on the HOMO level, since the inclusion of the LUMO led 
to severe sign-problems in QMC calculations. However, 
in the following calculations, the negligence of Coulomb 
interaction on the LUMO becomes a reasonable approx- 
imation since we choose the level spacing A much larger 
than the charging energy ( A » ^ ) such that the LUMO 
level is mostly empty. (See FIG. [T]) 

We model the el-ph coupling via the Jahn- Teller 



phonons of ajj 



if. 



ph,HL 



1,2) as 



_ (7) 
with the 2 X 2 JT coupling matrix [23'| given as F^^^ = 
.9ep<5'zj = QepO'x, with the Pauli matrices ct^ and 
ax- The second phonon (m = 2) makes direct elec- 
tronic transitions between the LUMO and HOMO with- 
out tunneling through reservoirs. The tunneling part is 
given as Ht.HL = -tY^aka Y^'ML^aka + h.c). The aver- 
age sign in the QMC calculations has been moderate at 
Sign = 0.5 -1.0. 

The conductance in FIG. [3] shows clear phonon excita- 
tion peaks at $ = Uph . The HOMO-LUMO level spacing 
is set at A = 15, much greater than any other energy 
scales. With the HOMO level at e<j = -1.1, the HOMO 
occupation number uhomo 0-57 at $ = for (a). For 
(c), Cd = —1.5 and uhomo ~ 0.48. The main features of 
the conductance are the ZBA and the peak at $ « ujph. 

To understand how the co-tunneling via phonon ex- 
citation arises, we study the QD spectral function as 
<i> — > ujph. In FIG. EJb), spectral functions are plotted 
for <I> = 0, ■ • • , 1 with the interval of 0.2. Curves are 
off-set for clarity. Destruction of the ZBA resonance is 

il,[ii|. 



At $ = 0, 



CT,j=l,2 



similar to the pure Anderson models 
the other dominant peak is the charge excitation peak, 
marked by A in (b). As grows, the peak A quickly mi- 
grates to the peak B at lu = ujph for phonon excitation. 
The peak B is not directly responsible for the conduc- 
tance peak at $ = ujph since the peak B is outside the 
transport energy window [— $/2,<I>/2] in Eq. ([5|). 

The co-tunneling transport is carried by a new emerg- 
ing resonance as indicated by peak C in FIG. El^b). As 
the ZBA peak disappears, another resonance peak in- 
side the transport energy window becomes stronger. At 
= ojphi the mismatch of electron Fermi energies from 
the source and drain is compensated by an emission of 
a phonon quantum. Then effectively the same electronic 
chemical potentials on both reservoirs seem to result in 
a Kondo-like phonon-assisted many-body resonance. 

In a pure phonon model in FIG. [D^c-d), the conduc- 
tance behavior remains qualitatively the same, but it 
showed a strong negative-differential-resistance (NDR) 
behavior near $/2 « tOph in the inelastic sequential tun- 
neling regime. As shown in a dashed line at $ = 2ujpin 
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FIG. 3: Differential conductance in the HOMO-LUMO model 
with Anderson- (Jahn- Teller) Holstein model, (a) Conduc- 
tance peaks for the Kondo anomaly and inelastic co-tunneling 
at $ = ujph- (b) Spectral functions for bias <I> = 0, • ■ ■ ,1 with 
the interval 0.2. Curves are shifted for clarity. With increas- 
ing <I>, the Kondo peak disappears and the Coulomb peak 
(peak A) shifts to phonon excitation energy (peak B). As 
$ — > (jjph, new resonance (peak C) emerges inside the volt- 
age window and contributes to the co-tunneling conductance 
peak at $ « u}ph- Dashed line is for f3 — 20. (c) Conductance 
in pure el-ph limit. A strong negative-differential-resistance 
(NDR) effect appears at <l?/2 = ujph in the sequential tunnel- 
ing regime, (d) Disappearance of the resonance (dashed line 
at <E>/2 = LUph) leads to the NDR. 



the spectral weight shifts to high frequency and the res- 
onance peak C is destroyed, which leads to the NDR. 
Similar polaronic effects to NDR in molecular systems 
have been reported previously [2J|. Although this be- 
havior is robustly reproduced at different parameters, it 
should be mentioned that this is the regime where the 
ansatz, Eq. ([U, starts to deviate from the QMC data 
significantly and vertex corrections may be necessary. 

Finally, nonequilibrium-induced multi-phonon modes 
are shown in FIG. 21 The phonon spectral functions are 
calculated by the same ansatz, Eq. ([4]), but with even 
integers 7. At $ = 0, the phonon frequency is highly 
renormalized from cuph, but phonons are mostly in the 
lowest quantum state. As $ increases, the spectral weight 
transfers to multiple phonon modes as previously pre- 
dicted ^E^. 

I thank useful discussions with F. Anders, H. van der 
Zant, R. Heary. I acknowledge support from the Na- 
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FIG. 4: (a) Phonon spectral function of HOMO/LUMO 
model for parameters of FIG. [3ja). At $ = 0, phonon is 
strongly renormalized. With increasing higher phonon 
quanta are created by absorbing the chemical potential dif- 
ference in the electron reservoirs, (b) Similar results for 
FIG.Hb). 
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